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Abstract. We present an optimization procedure for neutrino beams which could be produced at CERN 
and aimed to a set of seven possible underground sites in Europe with distances ranging from 130 km 
to 2300 km. Studies on the feasibility of a next generation very massive neutrino observatory have been 
performed for these sites in the context of the first phase of the LAGUNA design study. We consider specific 
scenarios for the proton driver (a high power proton driver at 4.5 GeV lor the shortest baseline and a 50 
GeV machine for longer baselines) and the far detector (a Water Cherenkov for the shortest baseline and a 
LAr TPC for longer baselines). The flux simulation profits of a full GEANT4 simulation. The optimization 
has been performed before the recent results on v e appearance by reactor and accelerator experiments and 
hence it is based on the maximization of the sensitivity on sin 2 2#i3 . Nevertheless the optimized fluxes have 
been widely used since their publication on the internet (2010). This work is therefore mainly intended 
as a documentation of the adopted method and at the same time as an intermediate step towards future 
studies which will put the emphasis on the performances of beams for the study of 8cp ■ 



PACS. 14.60.Pq Neutrino mass and mixing 



1 Introduction 

The feasibility of a European next-generation very massive 
neutrino observatory in seven potential candidate sites lo- 
cated at distances from CERN ranging from 130 km to 
2300 km, has been explored within the LACUNAE] design 
study [1 . In order of increasing distance from Geneva the 
considered sites have been Frejus (France) at 130 km, Can- 
franc (Spain) at 630 km, Caso (Italy) at 665 km, Sieroz- 
sowice (Poland) at 950 km, Boulby (United Kingdom) at 
1050 km, Slanic (Romania) at 1570 km and Pyhasalmi 
(Finland) at 2300 km. 

When coupled to very intense neutrino beams from 
CERN, large detectors hosted in such an underground site, 
could measure with high precision the mixing angle #13, 
and eventually determine the neutrino mass hierarchy and 
the existence of CP violation in the leptonic sector. 

The oscillation probability of the — > v e channel is 
shown as a function of the neutrino energy in Fig. [1] for 
the considered baselines. The energy of the first oscillation 
maximum spans a wide range of energies for the consid- 
ered baselines ranging from 0.26 MeV at 130 km up to 
4.65 GeV at 2300 km, the full sequence being {0.26, 1.27, 
1.34, 1.92, 2.12, 3.18, 4.65} GeV. 

This parameter is crucial for the optimization of the 
energy spectrum of the neutrino beam as it will be shown 
later. Neutrino spectra should cover the region where the 
oscillation effect is maximal with high statistics and low 
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intrinsic contamination of v e . The study of CP-violation 
requires to measure the oscillation probability as a func- 
tion of the neutrino energy, or alternatively to compare 
large samples of v e and v e CC events, and suffers in general 
from neutrino oscillation parameters degeneracies. The 
possibility to have a broad beam covering the second os- 
cillation maximum at lower energy is beneficial since it 
provides additional input useful to constraint the effects 
of mass hierarchy and the 5c p phase [2] and limits the 
impact of systematic errors on flux normalization by pro- 
viding spectral information. 

In this work we investigate two options for the pro- 
ton driver: a high power superconducting proton linac at 
4.5 GeV and a high power synchrotron at 50 GeV. Con- 
cerning the detector technology we consider a 440 kt Wa- 
ter Cherenkov for the 130 km baseline and the low energy 
proton driver and a 100 kt LAr Time Projection Cham- 
ber (LArTPC) at longer baselines with the high energy 
accelerator. Realistic designs have been proposed for these 
two detectors: the MEMPHYS 3^ and the GLACIER @] 
concepts. Previous studies on a high-energy super-beam 

5 and a low energy super-beam [7], [5], [5] [5] are avail- 
able. The simulation of fluxes is based on the GEANT4 

10 libraries and the optimization is performed separately 
for each one of the considered baselines. The guiding line 
of the optimization is the final achievable sensitivity on 
sin 2 2#i3. Furthermore a direct comparison of a high-energy 
and low-energy super-beam based on different accelerator 
scenarios has been done using of a coherent set of simula- 
tion tools. 
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Fig. 1. The neutrino oscillation probability V(y^ — > v E ) for the 
LAGUNA baselines. We assume the inverted hierarchy with 
sin 2 2013=0.05, Sep = and the presence of matter effects. 

2 Proton drivers and detectors 

Neutrino rates in conventional beams are at first approxi- 
mation proportional to the incident primary proton beam 
power, hence intense neutrino beams can be obtained by 
trading proton beam intensity with proton energy. So two 
basic approaches may be considered: a relatively low pro- 
ton energy accompanied by high proton intensity or a 
higher proton energy with lower beam current. 

In low energy neutrino beams the bulk of v e contam- 
ination comes from the it — > /i decay chain and only 
marginally from kaon decays (~ 10% of the v e at 4.5 GeV 
proton energy). This source of background can therefore 
be more easily constrained due to the correlation with the 
dominant flux component from direct pion decays. In 
addition, the horn and decay tunnel can be kept at rela- 
tively small scales. Finally in the sub-GeV region most of 
neutrino interactions are quasi-elastic. This final state al- 
lows an easy configuration for the calculation of the parent 
neutrino energy and can be cleanly reconstructed also in a 
water Cherenkov detector. The 7r° rejection from neutral 
current events also benefits from the low energy regime 
since photons are less collinear and energetic allowing the 
Cherenkov ring patterns generated by e + e _ to be more 
resolvable. In order to fulfill the condition of being on the 
first maximum of oscillation the baseline has to be con- 
formingly small and this offers the advantage of having 
a small suppression of the flux (~ L~ 2 ). Furthermore the 
determination of CP violation at small baselines is cleaner 
since there is almost no interplay with CP violating effects 
related to matter effects. 

On the other hand high energy super beams associ- 
ated to large baselines offer the possibility to study the 
neutrino mass hierarchy via the study of matter effects in 
the earth. The neutrino cross section scales about linearly 
with the energy allowing comparatively larger interaction 
rates at fixed flux. At high energy neutrino cross sections 
are free from the large theoretical uncertainties present 
in the low energy regime (nuclear effects, Fermi motion) 



which make the use of near detector compulsory for low 
energy super-beams. These effects also spoil the neutrino 
energy resolution. Furthermore the pion focusing is more 
efficient at high energy of the incident protons due to the 
more favourable Lorentz boost. We note finally that the 
chance to measure both first and second maxima increases 
with the baseline since, in general, the second maximum 
tends to fall at low energy where resolution and efficiency 
degrade. 

In the following we describe the assumptions for the 
proton drivers: 

— A conceptual design report (CDR2) exists for the high 
power super conducting linac (HP-SPL) [13]. It is fore- 
seen as a 4 MW machine working at 5 GeV proton 
kinetic energy. At a first stage it would feed protons 
to a fixed target experiment to produce an intense low 
energy (~ 400 MeV) super-beam. On a long time scale 
this machine could also be used to provide protons for 
the muon production in the context of a Neutrino Fac- 
tory. 

— The scenario for a high-power 50 GeV synchrotron 
(HP-PS2) has been initially proposed and discussed 
in [5]. A factor four in intensity is assumed compared 
to the baseline parameters defined by the PS2 working 
group [12] (1.2 • 10 14 protons with a cycle of 2.4 s). A 
value of 3 ■ 10 21 protons on target (p.o.t.) per year 
could be achieved by doubling both the proton inten- 
sity and the repetition rate. 



Parameter 


HP-SPL 


HP-PS2 


p kin. energy (GeV) 


5 


50 


repetition frequency (Hz) 


50 


0.83 


p per pulse (10 14 ) 


1.12 


2.5 


average power in 10 's (MW) 


4 


2.4 


p.o.t/year (lO^ 1 ) 


56 


3 



Table 1. Parameters of considered proton drivers. 



A summary of the assumed accelerator parameters is 
given in Tab. [1] For the far detector we concentrated on 
two designs: 

— The MEMPHYS water Cherenkov detector is envis- 
aged as consisting of 3 separate tanks of 65 m in di- 
ameter and 65 m height each (440 kt). Such dimen- 
sions meet the requirements of light attenuation length 
in (pure) water and hydrostatic pressure on the bot- 
tom PMTs. A detector coverage of 30% can be ob- 
tained with about 81.000 PMT of 30 cm diameter 
per tank. Based on the extensive experience of Super- 
Kamiokande, this technology is best suited for single 
Cherenkov ring events typically occurring at energies 
below 1 GeV. 

— GLACIER is a scalable concept for single volume very 
large LAr TPC with a mass of 100 kt. The powerful 
imaging will allow to reconstruct with high efficiency 
electron events in the GeV range and above, while con- 
siderably suppressing the neutral current background 
mostly consisting of misidentified 7r°s. 
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3 Optimization of the focusing system 

The optimization of the neutrino fluxes for the CERN- 
Frejus baseline with a Cherenkov detector and a 4.5 GeV 
proton driver has been studied extensively in [llj so in 
the following we will take the optimized fluxes obtained 
in that work and focus on the optimization of the focusing 
system for longer baselines assuming a LAr far detector 
and a 50 GeV proton driver. 

The focusing system is based on a pair of parabolic 
horns which we will denote as horn (upstream) and reflec- 
tor (downstream) according to the current terminology. 
This schema is the same which is being used for the NuMI 
beam. The target is modelled as a 1 m long cylinder of 
graphite (p = 1.85 g/cm 3 ) and a radius of 2 mm. Primary 
interaction in the target were simulated with GEANT4 
QGSP hadronic package. 

The optimization relies on a parametric model of the 
horn and reflector profiles. The horn radius as a function 
of the coordinate along the proton beam, r[z), has been 
parametrized as shown in the first row of Tabj2] in the 
three z domains [0, Z\\, [z\, Z2], [22, 23]. The model contains 
eleven shape parameters for each magnetic lens (a, b, c, d, 
a' , b', d , r, Zi, Z2, Z3) which reduce to nine after requiring 
continuity at the points Z\ and z-i (a, b, c, d, c', Zi, Z2, Z3). 
The layout of a typical configuration is shown in Fig. [5] 




z(m) z(m) 



Fig. 2. Layout of the horn/reflector. 



meson tracks per configuration was used. Fluxes were cal- 
culated with 20 energy bins from to 10 GeV. The sta- 
tistical fluctuations introduced by the size of the sample 
have been estimated by repeating the simulation for the 
same configuration several times with independent ini- 
tialization of the GEANT4 random number engine. The 
spread is enhanced by the presence of single events which 
can be assigned large weights. The spread on the param- 
eters Xl is of the order of 3-4%. The sensitivity limit was 
calculated with the GLoBES software [IB] fixing a null 
value for #13 and fitting the simulated data with finite 
values of sin 2 2#i 3 and 5c p sampled in a grid of 10 x 200 
points in the (Sep, sin 2 26*13) plane for Sep £ [0, 2vr] and 
sin 2 26»i 3 G [10" 2 ,10~ 4 ]. The 99% C.L. limit was set at 
the values corresponding to a A\ 2 of 9.21 (2 d.o.f.). The 
normal hierarchy was assumed in the calculation. We as- 
sumes running periods of 2 years in v mode and 8 years 
in v mode. The detector response is described in GLoBES 
by assigning values for the energy resolution, efficiency 
and defining the considered channels. The parametriza- 
tion of the MEMPHYS detector is the same which was 
used in [5]. The event selection and particle identification 
are the Super-Kamiokande algorithms results. Migration 
matrices for the neutrino energy reconstruction are used to 
properly handle Fermi motion smearing and the non-QE 
event contamination. The reconstructed energy is divided 
into 100 McV bins while the true neutrino energy in 40 
MeV bins from to 1.6 GeV. Four migration matrices for 
v e , f M , u e and are applied to signal events as well as 
backgrounds. The considered backgrounds are inter- 
actions misidentified as c , neutral current events and 
v R + D e intrinsic components of the beam. 

In the simulation of the GLACIER detector the consid- 
ered backgrounds are the intrinsic v e and v e components 
in the beam. Reconstructed neutrino energy was divided 
in 100 MeV bins from to 10 GeV. A constant energy 
resolution of 1 % is assumed for the signal and the back- 
ground. 

We followed two strategies in the optimization proce- 
dure which we will describe in the following subsections. 



In addition to the parameters related to the shape of 
the horn and the reflector, additional degrees of freedom 
are: the distance between the horn and reflector (Ahr), 
the length and radius of the decay tunnel (L tun , rt un ), the 
longitudinal position of the target (z ta r) and the currents 
circulating in the horn and the reflector (ih, ip.)- 

Following the approach already used in [11] for the 
optimization of the SPL-Frejus Super Beam, we introduce, 
as a figure of merit of the focusing, a quantity A defined as 
the 6cp-averaged 99 % C.L. sensitivity limit on sin 2 (26*13) 
(:=Xqq(Scp)) in 10~ 3 units 

1Q 3 /-27T 

A = — / \ 99 {S C p)d5cp (1) 
2tt Jo 

In the following we will denote the quantity A evaluated 
for a specific baseline L as Xl- A sample of 10 5 secondary 



r{z) 


y/ b c 


d 


/ z — a' 1 

y — ~ c 


z range 


[0, 21] 


[Z1,Z2] 


[Z2,Z S ] 


Par. 


horn 


refl. 


Par. 


horn 


refl. 


a 


85.7 


100 


d 


0.9 


3.9 


b 


7.0 


0.135 


r 


15 


40 


c 


0.2 


0.3 


Zl 


80 


97.6 


1 

a 


82.2 


100. 


Z2 


83.0 


104.8 


b' 


2.18 


0.272 




300 


300 


d 


0.2 


0.3 









Table 2. Analytic parametrization of the horn/reflector radial 
profile r(z) and central values of the parameters expressed in 
cm. r is the conductor outer radius. 
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Fig. 3. Dependence of A630 on zt ar f° r the fixed horn search. 



3.1 Fixed horn search 



column of Tab. EQ The 
and X' min shows that the tunnel 
optimization is particularly effective for the short base- 
lines for which the initial geometry was not appropriate. 
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In a first step we decided to fix the horn and reflector 
shapes (central values of Tab. [5]), the tunnel geometry 
(L tun = 300 m, r t un = 1-5 m) and the circulating cur- 
rents (200 kA). We then varied the relative positions of the 
horn, the reflector and the target. We define the distance 
between the center of the target and the most upstream 
point of the horn as z tar while we indicate with Ahr the 
horn-reflector distance. After having chosen the best point 
in this space we did a similar exercise in the decay tunnel 
parameter space (Ltun, rtun)- At first order these two cou- 
ples of parameters are expected to be weakly correlated so 
that doing the optimization in one pair of variables after 
fixing a specific choice for the other pair should not have 
a big impact on the final result. 

The variables (Ahr., z tar ) were sampled uniformly in 
the intervals [0, 300] m and [—1.5, 2.5] m respectively. Op- 
timal values were then chosen for each baseline. In general 
a marked dependence of A on the longitudinal position of 
the target (z tar ) is observed while variations of Ahr have 
a reduced impact. In Fig [3] we show, taking the baseline 
of 630 km as an example, the dependence of A630 on z tar 
after marginalizing on Ahr- For the 630 km baseline the 
optimal z tar lies around +0.5 m, while for Ahr. a value 
of 50 m was chosen. At this stage of the optimization the 
best values for A630 cluster between 1.4-1.5. The first two 
columns of Tab. 13. II give the Ahr and z tar pairs providing 
the best limit for each baseline (\ m i n , 3 rd column). 

After having fixed Ahr, z tar to the optimal values of 
Tab. l3.ll the tunnel length L tU n, previously fixed at 300 m, 
was sampled uniformly between [10, 500] m keeping r tun 
fixed at 1.5 m. The optimized values for L tun are given in 
Tab. ED (4*^ column). In the case of L = 630 km a gain of 
order 20% is visible in Fig. @] (left) decreasing L tun from 
a 300 to 75 m. 

The rtun was then sampled in [0,3] m having fixed the 
optimal tunnel length. The right plot of Fig. 0] shows that 
an improvement is obtained increasing r tun to 2 m. The 
optimized values for rt un are shown in Tab. O (5 th col- 
umn) . The values of A obtained after the tunnel optimiza- 
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Fig. 4. Dependence of A630 on L tU n and r tU n in the fixed horn 
search. 



3.2 General search 

In order to understand how limiting is the choice of fixing 
the shape of the conductors and the circulating currents 
we performed a high statistics general scan of the config- 
urations allowing also these parameters to vary. Also the 
parameters which were previously optimized were varied 
since in general changing the shape of the conductors and 
the currents we do not expect the same optimization to 
be valid any more. 

The shape parameters of the horn and the reflector 
were sampled with uniform distribution within 50 % of 
their central values given in Tab. [5] The other parameters 
were sampled with uniform distribution in the ranges of 
Tab. 01 The inclusive distributions of the input parame- 
ters were compared to the corresponding ones for the sub- 
sample with A <1.5 in order to pin down the variables 
which are effective in producing good results. Despite the 
smearing effect introduced by the simultaneous variation 
of many correlated variables, a visible trend is observed 
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L (km) 


Z tar ( m ) 


aTr m 




L tun ( m ) 


r tun ( m ) 


y 


630 


0.5 


50 


1.4-1.5 


75 


2 


1.3 


665 


0.45 


55 


1.4-1.5 


90 


2.2 


1.3 


950 





75 


1.3 


110 


2 


1.2 


1050 


-0.25 


4 


1.3 


200 


1 


1.3 


1570 


-0.3 


4 


1.2-1.3 


280 


1 


1.2 


2300 


-0.8 


4 


1.7 


400 


1.5 


1.6 



Table 3. Fixed horn shape search. Optimal values for Ahr, z tar , L tU n and r t 



Parameter 


interval 


Parameter 


Interval 


Ltun 


[200,1000] m 


itar 


2 mm 


'tun 


[0.8,2] m 


A H r 


[4,300] m 


Ztar 


[-2.5, 1.5] m 


tH, iR 


[150,300] kA 



Table 4. Focusing system parameters not related to the horn- 
reflector shapes. 




for ztar which exhibits a strong correlation with A. It is 
clear that putting the target more and more upstream 
with respect to the horn, is mandatory to get good exclu- 
sion limits, as far as the baseline increases. 

The correlation between the longitudinal position of 
the target with respect to the horn and the mean energy 
of the spectrum {{E„ )) is shown in Fig. [5] Putting 
the target upstream, high energy pions, which are typi- 
cally produced at small angles, are preferentially focused 
resulting in a high energy neutrino spectrum. 
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Fig. 5. Correlation between the longitudinal position of the 
target with respect to the horn and (E v ). 

The correlation between (E„ ) and A is shown in Fig. [6j 
In general the optimal energies tend to roughly follow the 
position of the first oscillation maximum (red vertical lines 
in the plots). Mean energies below 2 GeV are difficult to 
get with a 50 GeV proton beam. A possible solution, which 
has not be considered in this work, could be to go towards 
an off axis beam for baselines lower than 600 km. The 
horizontal blue lines show the lowest values for A obtained 
with the previous fixed horn shape search. 

The achieved performance is not drastically improved 
by the general search though some gain appears for L > 



v 




Fig. 6. Correlation between the figure of merit A and (E v ) 
(positive focusing). 

1000 km. Blue markers highlight the configuration pro- 
viding the best limit for each baseline. It turns out that 
the same configuration provides the best limit both for 
630 and 665 km and the same happens for 950-1050 and 
1570 km. Given the limited improvement, we decided to 
stick with the best candidates obtained with the fixed horn 
shape search. This choice is also motivated by the fact that 
choosing the configurations with the minimum A has the 
disadvantage of being sensitive to statistical fluctuations. 

4 Optimized fluxes 

The fluxes obtained with the optimized focusing setups 
according to the fixed shape search are shown in Fig. [7] at 
a reference distance of 100 km. Fluxes are publicly avail- 
able on the internet [T3]. The flux increase as the mean 
energy increases can be intuitively explained considering 
that high energy pions are easier to focus since they natu- 
rally tend to emerge from the target in the forward direc- 
tion and the neutrinos they produce have a higher chance 
to be in the far detector solid angle also thanks to the 
effect of the Lorentz boost. Un-oscillated interaction rates 
are given in Tab. [SJ It can be noted that considerable 
samples of r events becomes collectable with the fluxes 
optimized for the longer baselines. 

5 Conclusions 

As it was shown in [B], using the fluxes optimized with 
the procedure described above, the "discovery potential" 
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(km) 
130 
630 
665 
950 
1050 
1570 
2300 



41316 
36844 
38815 
37844 
51787 
26785 
17257 



174 
486 
516 
349 
314 
174 
110 



/ 

28 
28 
40 
148 
170 
377 



(%) 



0.42 
1.5 
1.5 
1.0 
0.64 
0.67 
0.67 



5915 
13652 
14287 
14700 
21728 
11184 
7577 



v run 








15 


1 


0.42 


157 


11 


2.0 


158 


11 


2.0 


107 


15 


1.3 


88 


65 


0.60 


47 


73 


0.57 


32 


172 


0.60 



Table 5. Charged current event rates with the optimized fluxes. Fractions are expressed as a percentage. Numbers are normalized 
to a detector mass of 100 kt and a running time of one year corresponding to 3 • 10 21 p.o.t. for the 50 GeV proton driver and 
56 • 10 21 p.o.t. for the 4.5 GeV option. 




Fig. 7. Neutrino fluxes at 100 km for the systems optimized 
with the fixed horn shape search. The energies of the oscillation 
maximum for each baseline (Fig. [TJ are indicated with vertical 
lines having the same color as the corresponding spectrum. The 
integral fluxes in units of 10 15 t' M /100 m 2 /year are 0.38, 1.59, 
1.81, 2.69, 3.56, 3.93 and 4.48 in order of increasing baseline. 

for 0i3 turned out to be, at first order, almost indepen- 
dent of the baseline. Performances of high- and low-energy 
super-beams are comparable if we assume for both a 5% 
systematic error on the fluxes. Concerning the high-energy 
super beam, better results are obtained for intermediate 
baselines from 950 to 1570 km even though the difference 
is not marked. This merit factor, despite being obsolete 
after the recent experimental results, is probably still a 
reasonable indicator of the precision with which #13 could 
be measured with these configurations. 

This result could be achieved by a systematic tuning of 
a few basic parameters of the focusing system: the horn- 
reflector distance, the target position and decay tunnel 
geometry. 

An exhaustive discussion of the physics potential in 
terms of CP violation and mass hierarchy obtainable with 
the fluxes whose optimization is described in this work, 
has been recently developed in [16], [17], [18], [19] and 

By adopting a suitable re-definition of the figure of 
merit, the approach followed in this study could in the 
future be specialised to the need for optimal sensitivity 



on the CP violating effects under the light of the recent 
measurement of #13. 
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